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ABSTRACT 


Crocco's three-dimensional nozzle admittance theory is extended to 
he applicable when the amplitudes of the combustor and nozzle oscillations 
increase or decrease with time. An analytical procedure and a computer 
program for determining nozzle admittance values from the extended theory 
are presented and used to compute the admittances of a family of liquid- 
propellant rocket nozzles . The calculated results indicate that the nozzle 
geometry, entrance Mach number and temporal decay coefficient significantly 
affect the nozzle admittance values. The theoretical predictions are shown 
to be in good agreement with available experimental data. 
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INTRODUCTION 


The interaction between the pressure oscillations inside an unstable 
rocket combustion chamber and the wave motion in the convergent section of 
the exhaust nozzle can have a significant effect on the stability charac- 
teristics of the rocket motor and is an important consideration in analytical 
studies concerned with the prediction of the stability of liquid-propellant 
rocket engines. This report is concerned with the investigation of this 
interaction. 

To determine the stability of a liquid-propellant rocket engine, the 
equations describing the behavior of the oscillatory flow field throughout 
the rocket motor must be solved. To simplify the problem, it is convenient 
to analyze the oscillations in the combustion chamber and the nozzle separately. 
For such an analysis, the combustion chamber extends from the injector face to 
the nozzle entrance as shown in Fig. 1. All the combustion is assumed to take 
place in the combustion chamber where the mean flow Mach number is generally 
assumed to be low. On the other hand, no combustion is assumed to take place 
in the nozzle and its mean flow Mach number increases from a low value at the 
nozzle entrance to unity at the throat. Downstream of the throat the flow 
is supersonic and disturbances in this region cannot propagate upstream and 
affect the chamber conditions. Therefore, in combustion instability studies 
it is only necessary to consider the behavior of the oscillations in the 
converging section of the nozzle since only these oscillations can influence 

the conditions in the combustion chamber. 

1 2 

The nozzle admittance 5 is the boundary condition that must be 
satisfied by the combustor flow oscillations at the nozzle entrance. Defined 
as the ratio of the axial velocity perturbation to the pressure perturbation 
at the nozzle entrance, the nozzle admittance can also be used to determine 
whether wave motion in the nozzle under consideration adds or removes energy 
from the combustor oscillations. Furthermore, this boundary condition 
influences the structures and resonant frequencies of the natural modes of 
the combustor under investigation. 

To theoretically determine the nozzle admittance, the equations 
which describe the behavior of the waves in the convergent section of the 
exhaust nozzle must be solved. These equations have been developed by 
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Crocco and were solved numerically to obtain admittance values for one- 
and three-dimensional oscillations. These values were tabulated over a 
wide range of frequencies and entrance Mach numbers for a specific nozzle 
geometry. By applying the scaling technique developed in Ref. 2, the 
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admittances of related nozzles can be determined. It was pointed out, 
however, that interpolation of the tabulated values can result in large 
errors in the predicted nozzle admittances; furthermore, the accuracy of 
the scaling procedure is open to question. In addition, Crocco’ s theory 
is only applicable to constant amplitude periodic wave motions, and in its 
present form it cannot be applied to cases where the amplitude of the 
oscillations varies in time. 

In this report, the equations needed for computing the nozzle admit- 
tance are presented and their solutions are outlined. Crocco' s theory is 
extended to account for wave-amplitude variation with time. Typical 
theoretical, predictions are shown and compared with available experimental 
data. The effects of the nozzle geometry and chamber Mach number on the 
nozzle admittance are presented in plots showing frequency dependence of the 
real and imaginary parts of the nozzle admittance. The effects of the decay 
coefficient are also assessed. A manual describing the use of the computer 
program which calculates nozzle admittance values along with a program 
listing is presented in the appendix. 

SYMBOIS 


A, B, C 
c 

% 

i 

J m 

K(M>t) 


variable coefficients defined below Eq. (14) 
nondimens ional speed of sound, c*/c* 
unit vectors 

Bessel function of the first kind of order m 
a function having the following space and time dependence: 


M 


J [s (J ! -) 2 ]e 1 
mL mn\f / J 


icut ± im9 


' w 

Mach number at the nozzle entrance 
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m 

n 

P 

q 

r 

r 

cc 


S 

S mn 


t 

u 

v 

w 

y 


z 

Y 

C 

e 

9 

X 

P 

T 

<P 

I 


0) 


number of mode diametral nodal lines 
number of mode tangential nodal lines 
nondimensional pressure, p*/p* 
nondimens ional velocity, q*/c* 

nondimens ional radius, r*/r* 

nondimensional radius of curvature at the nozzle entrance, 

r */r* 
cc' c 

nondimens ional radius of curvature at the nozzle throat, 

r f/r* 
cv c 

nondimens ional frequency, cu*r*/5* 
the nth root of the equation 


dx 


nondimens ional time, t*c*/r* 
nondimensional axial velocity component, u*/c* 
nondimens ional radial velocity component, v*/c* 
nondimensional tangential velocity component, w*/c* 
irrotational specific nozzle admittance defined in Eq. (13) 

/ 

- - u * -- u 

y = p* c * — = Y PC — 

P '* v' 

nondimensional axial coordinate, z*/r* 
ratio of specific heats 

a function used to compute the nozzle admittance; defined below 
Eq. ( 13 ) 

tangential coordinate, radians 
nozzle half-angle, degrees 

nondimensional temporal decay coefficient, A*r*/c* 
nondimensional density, p*/p* 

a function used to compute the nozzle admittance; t = l/Q 
nondimensional steady state velocity potential, cp*/c*r* 
a function describing the cp- dependence- of the radial velocity 
perturbation 

— — 2 

nondimensional steady state stream function, yp^) q(cp) r 
nondimensional frequency, au*r*/c* 
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Subscripts : 


c evaluated at the chamber wall 

i imaginary part of a complex quantity 

o stagnation value 

r real part of a complex quantity 

th evaluated at the nozzle throat 

w evaluated at the nozzle wall 

-* vector quantity 


Superscripts: 

1 perturbation quantity 

— steady state value. 

* dimensional quantity 

ANALYSIS 


Derivation of the Wave Equations 

2 

The equations used by Crocco to compute the nozzle admittance will 
be developed from the conservation equations. To keep the problem mathe- 
matically tractable and yet physically meaningful, the following assumptions 
were employed. 

(1) The nozzle flow is a calorically perfect gas consisting of 
a single species. 

(2) Viscosity and heat conduction are negligible. 

(3) The steady state flow is one -dimensional; this assumption 
implies that the nozzle is slowly converging. 

(4) The amplitudes of the waves are small so that only linear 
terms in the perturbed quantities need to be retained in the 
conservation equations. 

(5) The oscillations are assumed to be irrotational. 

Using these assumptions, the equations of motion in nondimens ional 
form become 

Continuity 

+ V.(pq) = 0 (1) 
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Momentum 


Sq 

]_~2 1 . * 

^_ + ^Vq = - -Vp (2) 

2 / V 

and, from the isentropic conditions, c = p/p and p = p . 

To obtain the linearized wave equations, the dependent variables 
are expressed in the following form: 


q = q + q', p = p + p', p = P + p' 
— * 


( 3 ) 


Substituting these expressions into Eqs. (1) and (2), neglecting all non- 
linear terms involving primed quantities, and separating the resulting system 
of equations into a set of steady state equations and a set of unsteady 
equations yield the system of steady state equations : 


V* (pq) = 0 : 


-2 

c 


rY - 1 


= 1 


y - 1 


-2 

q ; 


-Y 

p = p 


( 4 ) 


and the following system of ■unsteady linear equations that describe the wave 
motion: 


^ + V.(qp' + pq') = 0 (5) 

(6) 

(T) 

To simplify the application of the boundary conditions at the nozzle 
walls, these wave equations are solved in the orthogonal coordinate system 
shown in Fig. 1 . In this coordinate system the steady state velocity 
potential cp replaces the axial coordinate z, the steady state stream function 
ijt replaces the radial coordinate r and the angle 6 is used to denote azimuthal- 
variations. Using this coordinate system the velocity vectors can be expressec 
as follows: 


-4 , /\ 

-gr- + VCq-q') = - v(£_J 

YP 


5 



5 = i(«P)e v 


$ = + v'^ + w' % 

Using the definitions of the steady state velocity potential and stream 
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function for a one -dimensional mean flow, it can he shown that 





= |p('P)i('P)r 2 


Rewriting Eqs . (5) and (6) in the (cp,i|r,9) coordinate system yields the 
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following system of equations : 


Continuity 



+ * 


a ( £l , ^ 

5cp \ p + q. / 



H 5 (rw') _ o 
2i|r 59 “ u 


Momentum 

cp-component 



ijr- component 



( 8 ) 


( 9 ) 


( 10 ) 


9 -component 



( 11 ) 


Equations (7) through (11) constitute a system of five equations in the five 
unknowns p / /p, u / /q, v'/rpq, rw / , and p^/yp* These equations are solved 
by the method of separation of variables and the solutions are 


6 



q 


— ~ = ■gd K (t? 9 jt) J 

rpq ^ 


rw' = *(«P) 3g- L K(i|f } e,t) 


[i(u> - iX)f (cp) + q 2 (cp) - J K(t,9 ,t) 


^[i(«J - U) f(c p) + q 2 (cp) ^l] K(t,0 ,t) 


K(t,e,t) = 


(j S (i)H cos mee 3 "^ 13 

I mL mn\\|r / J 


j iv(ff> ±ioS ^ ■ iX)t 


for standing waves 


for spinning waves 


These solutions identically satisfy the momentum and energy equations. 
Substituting these solutions into Eq. (8) and eliminating variables give 
the following differential equation for the function § : 

-2,-2 -2v d 2 § -2f l dq 2 / ..sldl 

q (c - q) — - q + 21 ^ - iX) J — • 

dcp c 


+ [(«J 


2-2 

_p _p q o 

. ^ \ 2 Y - 1 • f . ^ \ q dq b mn ~| 

e r 

w 


the formula' 


The function § can be related to the specific acoustic admittance by 

-n _2 


, — u YpcQ 

y = Ypc — = - - Z g — — c-a - — 

p q C + i(«> - iX) 
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*1 (][$ 

where C = j rg-. Using the definition of £ and Eq. (12) , the following 
differential equation for £ is derived: 


d£ B c + ,2 C 

dcp ” A b A 


( 14 ) 


where 


A 


- 2,-2 

q (c 



B 





a 2 - 2 
S c 

ran 


w 


- i(H> - iA) 



Equation (14) is a complex Riccati equation which must be solved 
numerically to obtain £. Once the value of £ is determined at the nozzle 
entrance, the nozzle admittance can be computed directly from Eq. (13) • 
Inspection of Eq. (14) shows that the value of £ depends upon its coefficients 
A, B, and C which in turn depend upon U), A, S , and the space dependence of 
q and c in the nozzle. The behavior of q and c in the nozzle can be computed 
once the value of Y and the nozzle contour are specified. 

To determine £ for given values of ou. A , and V and a specific 

nozzle contour, Eq. (14) must be integrated numerically. A major difficulty 

which can occur during this integration is that £ becomes unbounded whenever 

§ approaches zero, which causes numerical difficulties in the integration 
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scheme. Crocco and Sirignano noted that this phenomenon occurred for low 
Mach numbers and high values of At these Mach numbers and frequencies 

they developed asymptotic solutions for £. 

Instead of using the asymptotic solution, an exact numerical solution 
is obtained in this study. The problem is resolved by introducing a new 
dependent variable 


T 


1 l 
£ ~ d£ 
dcp 
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As § approaches zero and the magnitude of £ becomes large, T becomes smalls 
Introducing the definition of f into Eq. (14) gives the following Riccati 
equation for t 


It b 

dcp A 


T T 2 = 1 
A 


(15) 


At those regions where £ becomes unbounded, Eq. (15) is integrated instead of 
Eq. (14) . 


Method of Solution 

To obtain the nozzle admittance from Eq. ( 13 ) , values of £ and t are 
computed by numerically integrating Eq. (l4) or (15) . To evaluate the coeffi- 
cients A, B, and C, a differential equation that describes the variations of 
the steady state velocity in the subsonic portion of the nozzle must be derived. 
Differentiating the continuity equation 

5r2 « - 4h r th 4h = constant < 16 > 


- 2 - 2 

where = c, ^ = 2/ (y + 1) , and using Eq. (4) yield the following differential 
equation 


dq 2 = 1 

dr dr/dq 2 



(17) 


Using Eq. (17) and the specified nozzle contour in terms of r(z) , the quantity 
dq/dcp can be obtained from the relationship 

-2 2 - 

dq _ dq dr dz _ ^ JU dr 

dcp ~ dr dz dcp _ dr dz ^ ' 

_ p 

Once q is known the corresponding value of c (cp) can be obtained by 
use of Eq. (4) . To evaluate dr/dz in Eq. (.18) , the nozzle contour shown in 
Fig. 2 is used. Starting at the combustion chamber the contour is generated 
by a circular arc of radius r cc turned through an angle 9^, the nozzle 
half- angle. This arc connects smoothly to a straight line which is inclined 
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at an angle 9^ to the nozzle axis . This straight line then joins with another 
circular arc of radius r . which turns through an angle 9^ and ends at the 
throat* Using this nozzle contour, in regions I, II and III of Fig. 2 


dr 

dz 


I 


[ 2r ot (r 


th 


) - (r 


l 

\2-|2 

T th> ] 


ct 


+ r 


th 


- r 


dr 

dz 


dr 

dz 


= - tan 9. 


II 


III 


2-12 


[2r cc (l - r) - (1 - r) ] 

1 > r - r 
cc 


Utilizing the appropriate expression for dr/dz, Eq. (18) can now he solved 

simultaneously with Eq. (14) or (15) to determine the nozzle admittance. 

The numerical integration of these equations must start at some 

initial point where the initial conditions are known. Since the equation 
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for C is singular at the throat , the integration is initiated at a point 
that is located a short distance upstream of the throat. The needed initial 
conditions are obtained by expanding the dependent variables in a Taylor 
series about the throat. To obtain this Taylor series, its coefficients 


C(0) = £ 0 and m 


dC 


dcp 


must be evaluated at the throat where cp = 0. 

|cp * 0 

These coefficients are evaluated by substituting the series 

C = C o + C-jtp + . . . 

into Eq. (14) and taking the limit as cp -» 0. The results are 


Cq ~ C(0) - B 


0 


fe l ~ dcp 


cp = 0 


- W^) - - c J/< a x - v 


where 
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c o =G 


cp = 0 


[ (m -U) g -i g ( y - _ %(2 )] 


(y + l)/r,, F r. , 
v ' th ct th 


B~ = B 


0 cp =0 y + 1L 


4 T 1 ./ 

r 1 + x(to - i\) 


7r,, r 7 

th ct 


B 

B 1 dcp 


4 


6 + Y + i 2 (cju - iX) J 


cp = o Y + 1 3r th r ct /r th r ct 


A 

1 dcp 


4 


tp = ° (Y + D/r th r ct 


p _ ac 
°1 " dcp 


Y - St- 


mn 


. >Y + X/L 2/ 

cp = 0 r. , /r, , r , 

th th ct 


i(m - 11) 

3r th r ct 


(6 + y) 


The following relations are used in the evaluation of the above quantities: 


cp = 0 Y + 1 


dq 2 

dcp 


cp = 0 (y + l)/r,,r 


th ct 


Once Cq and £^ are known, the initial condition at cp = cp^ is obtained from 
the expression £(cp^) = £q + 

The numerical solution is obtained by use of a modified Adams 
predictor-corrector scheme, and employing a Runge-Kutte scheme of order 
four to start the numerical integration. Initially, Eqs. (14) and (18) 
are integrated to determine £; if the magnitude of £ exceeds a specified 
value at which numerical difficulties can occur, the integration of Eq. (14) 
is terminated. Using the value of £ at that point, t is computed and the 
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integration proceeds using Eq. (15) . Similarly, should the magnitude of x 
become excessively large, the integration of Eq. (15) is terminated, C is 
computed from the value of t at that point, and the integration proceeds 
using Eq. (14) . This process is repeated until the nozzle entrance is 
reached. A computer program utilizing this procedure has been •written in 
FORTUM V for use on the UNIVAC 1108 computer and it is presented in the 
Appendix. 


RESULTS AND DISCUSSION 

Using the previously mentioned computer program, theoretical values 

of the real and imaginary parts of the nozzle admittance have been computed 

for several nozzle configurations having contours similar to the one presented 

in Fig. 2. In these computations the radii of curvature, r^ c and r^ , are 

assumed to be equal. The admittance values are presented as functions of the 

nondimensional frequency S in Figs. 3 through 9 where they are compared with 

available experimental data obtained from Ref . 3* In these figures, the 

frequency has been nondimens ionali zed by the ratio of the steady state speed 

of sound at the nozzle entrance to the chamber radius r . 

c 

Admittances for Longitudinal Modes 

Longitudinal- type instabilities in general occur in the range of S 
from 0 to approximately 1.8 which is in the vicinity of the cutoff frequency 
of the first tangential modes. The cutoff frequency of a particular transverse 

mode is S V (1 - W) where S is the transverse mode eigenvalue and the sub- 
mn mn 

scripts m and n respectively denote the number of diametral nodal lines and 
the number of tangential nodal lines. Values of S mn are given in Table 1 for 
several values of m and n. 

For longitudinal modes good agreement exists between the experimental 
and theoretical values of the real and imaginary parts of the admittance as 
shown in Figs. 3 through 5. The effect of changing the nozzle half-angle is 
presented in Fig. 3 for a nozzle with an entrance Mach number M of 0..08 and 
r /r =s 0.44. The data indicate that increasing 0, increases the frequency 
at which the real and imaginary parts of the admittance attain maximum values. 
These data also indicate that the assumption of a one- dimensional mean flow 
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Table 1. Values of Transverse Mode Eigenvalues; S. 


Transverse Wave Pattern m 

Longitudinal 0 
First Tangential (IT) 1 
Second Tangential (2T) 2 
First Radial (1R) 0 
Third Tangential (3'T) 3 
Fourth Tangential (4 t) 4 
First Tangential, First Radial (1T,1R) 1 
Fifth Tangential (5T) 5 
Second Tangential, First Radial (2T, 1R) 2 
Second Radial (2R) 0 


mn 
n 

0 

0 

0 

1 

0 

0 

1 

0 

1 

2 


mn 

0 

1.8413 

3.0543 

3.8317 
4.2012 
5.3175 
5.3313 
6 . 4154 
6.7060 
7.0156 


used in the development of the theory appears to be valid. Even for nozzles 
with half-angles as high as 45 degrees, for which it has been shown that the 
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mean flow is two-dimensional? the experimental and theoretical nozzle admit- 
tance values are in good agreement. 

Examination of Fig. 4 shows that the entrance Mach number M has a 
significant effect on the admittance values for 9^ 05 15 degrees and r cc / r c = 
0.44. However, increasing the nozzle half-angle appears to decrease the 

influence of the entrance Mach number, and for 9 = .45 degrees variations in 

3 X 

M has little effect. The dependence of the nozzle admittance upon the radius 

of curvature for a nozzle with M = 0.16 and 9^ = 30 degrees is shown in Fig. 5. 

The data presented in Figs. 3 through 5 show that for longitudinal 

modes the real part of the nozzle admittance is always positive. As indicated 
1 2 

by Crocco 5 positive values of the real part of the nozzle admittance imply 
that the nozzle removes acoustic energy from the combustor wave system which 
implies that the nozzle exerts a stabilizing influence upon the chamber 
oscillations . 

In combustion instability analyses of liquid-propellant rocket motors, 
it is often assumed that the nozzle is short. This assumption implies that 
the nozzle length and throat diameter are much smaller than the chamber length 
and diameter so that the wave travel time in the nozzle is much shorter than 
the wave travel time in the chamber. For a short nozzle the real and imaginary 
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parts of the admittance are independent of frequency and are given by the 
expressions'* 



These theoretical short nozzle admittance results do not agree with the 
results obtained for typical liquid rocket nozzles presented in Figs. 3 
through 5. The disagreement is especially evident for nozzles with low 
values of 9^, which imply that the nozzle is long, and for high values of 
S where the wave length of the oscillation becomes of the same order of 
magnitude as a characteristic nozzle dimension. 

Admittances for Mixed First Tangential-Longitudinal Modes 

The mixed first tangential-longitudinal modes are those three- 
dimensional modes which exist between the cutoff frequencies of the first 
tangential (S ~ 1.8) and second tangential (S a 3.0) modes. Theoretical 
and experimental nozzle admittance data for these modes are presented in 
Figs. 6 through 8. 

In Fig. 6 the influence of the nozzle half- angle on the admittance 

values is shown. The theoretical and experimental results are in good 

agreement and they indicate that increasing 9^ increases the frequency at 

which the real and imaginary parts of the admittance reach maximum values. 

The effect of Mach number on the admittance values is presented in 

Fig. 7 for 9 n = 15 degrees and r /r ■ 0.44. Mach number effects are 

especially significant at the higher frequencies. However, as shown in 

Ref. 3, increasing the nozzle half-angle decreases the dependence of the 

admittance values on the Mach number. The effect of changing the radii of 

curvature on the admittance values is presented in Fig. 8. 

The results presented in Figs. 6 through 8 show that for mixed 

first tangential- longitudinal modes the real part of the nozzle admittance 

can be negative which means that the nozzle radiates wave energy back into 

the combustor; this process exerts a destabilizing influence on the oscil- 

2 

lations in the chamber. These negative values occur only for three- 

2 

dimensional modes and, as shown by Crocco, their cause can be traced to 

the term involving S in Eq. (12) . For longitudinal modes, for which S 

mn mn 
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is zero, the real part of the nozzle admittance is always positive, and for 
those modes the nozzle always exerts a stabilizing influence upon the combustor 
oscillations . 

Effect of Decay Coefficient upon Admittance Data 

The nozzle admittance theory has been modified to include the effects 
of a temporal decay coefficient, X. Typical results are shown in Figs. 9 
and 10 for values of X of -0 .05, 0, and 0,05. These results indicate that 
varying X affects both the real and imaginary parts of the admittance. There- 
fore, the decay coefficient should be included in the nozzle admittance 
computations when the oscillations are not neutrally stable. 

SUMMARY AID CONCLUSIONS 


The equations necessary to determine the nozzle admittance for one- 
and three-dimensional oscillations have 'been developed. The analytical 
approach used in solving the nozzle wave equations is outlined and employed 
to obtain nozzle admittance data for typical nozzle configurations. These 
data show the dependence of the nozzle admittance values upon nozzle geometry, 
nozzle Mach number, mode of oscillation, and the temporal damping coefficient. 

The results can be summarized as follows for longitudinal and mixed 
first tangential- longitudinal modes . Decreasing the nozzle length by increas- 
ing the nozzle half-angle and Mach number or by decreasing the throat and 
entrance radii of curvature decreases the frequency dependence of the nozzle 
admittance. Good agreement exists between the theoretical predictions and 
available experimental data. However, the nozzle admittance values for typical 
liquid rocket nozzles are not in agreement with the values obtained from short 
nozzle theory. Including the effects of a temporal damping coefficient in 
the nozzle admittance computations changes the admittance values. Therefore, 
when the oscillations are not neutrally stable, the temporal decay coefficient 
should be accounted for in the computations . 
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Figure 1. Typical Mathematical Model of a Liquid Rocket Engine 



Figure 2 , Nozzle Contour 
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Figure 3 • The Effect of Nozzle Half -Angle on the Theoretical and 

Experimental Nozzle Admittance Values for Longitudinal Modes 
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Figure 4 . The Effect of Entrance Mach Humber on the Theoretical and 

Experimental Nozzle Admittance Values for Longitudinal Modes 
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Figure 5» The Effect of the Radii of Curvature on the Theoretical and 
E:xperimental Nozzle Admittance Values for Longitudinal Modes 












Effect o 
irimental 
;ential-L< 






s 

Effect of the Radii of Curvature on the Theoreti 
irimental Nozzle Admittance Values for Mixed Firs 
;ential- Longitudinal Modes 









Figure 9. Effect of the Temporal Decay Coefficient on the 

Theoretical Nozzle Admittance Values for Longitudinal Modes 
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APPENDIX 

COMPUTER PROGRAM USED TO DETERMINE THE IRROTATIONAL NOZZLE ADMITTANCE 

The computer program for calculating the irrotational nozzle admit- 

2 

tance from Crocco's theory which is extended to account for temporal damp- 
ing is written in FORTRAN V interpretive language compatible with the 
UNIVAC 1108 machine language compiler. This program consists of seven 
routines - the main or control program and six subroutines. The names of 
the routines are listed in Table A-l in sequential order. The FORTRAN 
symbols used in these routines and their definitions are presented in Table 
A-2 in alphabetical order. The input parameters necessary for the admit- 
tance computations must be specified in the main program and are listed in 
Table A-3- The output parameters and their definitions are listed in Table 
A-k. A detailed flow chart of the computer program is shown in Fig. A-l, 
and the program listing and sample output are presented in Tables A-5 and 
A- 6, respectively. 

This computer program has been written to predict nozzle admittances 
for nozzle contours shown in Fig. 2. The run time required depends upon 
the number of admittance values desired and the nozzle length. To obtain 
10 admittance values at different frequencies for the nozzles investigated 
in this study, one to two minutes of run time on the UNIVAC 1108 computer 
are required . 
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Table A-l. List of Subroutines in the Computer Program Used 
to Determine the Irrotational Nozzle Admittance 


Subroutine 

Description 

MAIN 

Specifies the nozzle geometry and 
operating conditions in the converging 
section of the nozzle 

NOZADM 

Specifies initial conditions at the 
throat, computes the final nozzle admit- 
tance values, and contains all output 
formats 

RKTZ 

Uses the Runge-Kutta of order four to 
obtain initial values for the modified 
Adams integration routine 

RKZDIF 

Computes the differential element in 
the converging section of the nozzle 
used to solve Eq. (1*0 

KKTDIF 

Computes the differential element in the 
converging section of the nozzle used to 
solve Eq. (15) 

ZADAMS 

Numerically integrates Eq. ( 1*0 using 
the modified Adams numerical integration 
Scheme 

TADAMS 

Numerically integrates Eq. (15) using 
the modified Adams numerical integration 
scheme 
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Table A-2 . Definition of FORTRAN Variables 
(Page 1 of 4) 


Variable 


Definition 


A 

A(5) 

AF 

ANGLE 

AIR 

BI 

BR 

BOI 

BOR 

B1I 

B1R 

C 

Cl 

CM 

C0R(5) 

CR 

COI 

COR 

C1I 

C1R 

DP 

DP (5) 
DR 
DU 
DWG 


Real coefficient A of Eqs . (l4) and (15) 

Coefficients of the Runge-Kutta formulas of order four 
Nondimensional temporal damping coefficient \ 

Nozzle half -angle, degrees 

Derivative of the coefficient A evaluated at the throat 
Imaginary part of the coefficient B in Eqs . (l4) and ( 15) 
Real part of the coefficient B in Eqs. (14) and (15) 

Value of BI at the throat 

Value of BR at the throat 

Derivative of BI evaluated at the throat 

Derivative of BR evaluated at the throat 

2 

Nondimensional speed of sound squared, c 

Imaginary part of the coefficient C in Eqs. (l4) and (15) 

Mach number at the nozzle entrance 

Formula for the corrector in the modified Adams inte- 
gration routine 

Real part of the coefficient C in Eqs. ( 14) and (15) 

Value of Cl at the throat 

Value of CR at the throat 

Derivative of Cl evaluated at the throat 

Derivative of CR evaluated at the throat 

Integration stepsize 

Derivative used in the corrector formula in the modified 
Adams integration routine 

Derivative of the local wall radius with respect to 
axial distance 

_2 

Derivative of the nondimensional velocity q with respect 
to the wall radius r 

Increment of the nondimensional frequency 00 
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Table A-2. Definition of FORTRAN Variables 
(Page 2 of 4) 


Variable 


Definition 


dy(5,-4) 

F 

fz(4,5) 

FI 

F2 

GAM 

G(5) 

H 

I 


IQ 

IQZ 

J 

JOPT 

K 

N 

NU 

NWC 

P 

PARG 

PHII 

PHIR 


Derivative used in the modified Adams integration scheme 

Constant given as q/yp evaluated at the nozzle entrance 

Derivative used in the Runge-Kutta method 

Lumped parameter determined by the conditions at the 
throat 

Lumped parameter determined by the conditions at the 
throat 

Ratio of specific heats y 

Dependent variable in the Runge-Kutta integration routine 
Integration stepsize 
Integer counter 

Integer constant. If IP = 0 the nozzle admittance is 
output. If IP ^ 0 the amplitude and phase of the 
pressure oscillation are output along the length of the 
nozzle 

If IQ == 2, the integration of Eq. (15) for t is complete 

= 1: Eq. (15) lor t is integrated 

= 2: Eq. (14) for £ is integrated 

Integer variable 

= 1: Eq. (15) for T is integrated 

= 2: Eq. (l4) for £ is integrated 

Integer variable 

Integer variable 

Number of differential equations to be solved by the 
Runge-Kutta or the modified Adams integration routine 

Number of frequency points 

Value of the steady state velocity potential 
Phase of the pressure oscillation in the nozzle 
Imaginary part of $ 

Real part of $ 
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Table A-2. Definition of FORTRAN Variables 
(Page 3 of 4) 


Variable 


Definition 


PI 

PMAG 

PR 

PRED(5) 

Q 

QBAR 

R 

RCC 

RCT 

RHO 

RT 

Rl 

R2 

SRTR 

SVN 

SVNR 

SYI 

SYR 

T 

TDN 

TI 

TMAG 

TPI 


Imaginary part of the pressure oscillation 

Magnitude of the pressure oscillation 

Real part of the pressure oscillation 

Predictor formula for the modified Adams integration 
routine 

Y + 1 

2 % - 1 ) 

Constant given as (r^/4) (■ -— 

K ondimens i onal steady state velocity q 
Local wall radius r 

Ratio of the radius of curvature at the nozzle entrance 
to the radius at the nozzle entrance 


Ratio of the radius of curvature at the throat to the 
radius at the nozzle entrance 

N ondimens i onal , steady-state density p 

Nondimensional throat radius 

Nondimensional radius at the entrance to Section 2 of 
the converging portion of the nozzle 

Nondimensional radius at the entrance to Section 3 of 
the converging portion of the nozzle 

Constant give as / r ,, r / r 

tn cc c 

S mn 

S im r a/ r th 

Imaginary part of the specific admittance y 

Real part of the specific admittance y 

Nozzle half- angle, in radians 

Inverse of the square of the magnitude of Q 

Imaginary part of t 

Magnitude of t 

Derivative of TI with respect to cp 
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Table A- 2® Definition of FORTRAN Variables 
(Page 4 of 4) 


Variable 


Definition 


TPR 

TR 

TZ 

T2 

U 

UZ 

W 

WC 

X 

Y(5) 

YI 

YR 

ZDN 

ZI 

ZMA.G 

ZPI 

ZPR 

ZR 

ZOI 

ZOR 

Z1I 

ZlR 

Z2 


Derivative of TR 'with respect to cp 
Real part of t 

Value of cp at the nth integration point 
Square of the magnitude of t 

_2 

Steady state velocity squared, q 

Dependent variable in the Runge-Kutta integration scheme 

Nondimensional frequency S 

Nondimensional frequency ou 

Value of cp at the nth integration point 

Dependent variable used in the modified Adams integration 
scheme 

Imaginary part of the irrotational nozzle admittance 
defined by Crocco in Ref. 2 

Real part of the nozzle admittance defined by Crocco 
in Ref. 2 

Inverse of the square of the magnitude of £ 

Imaginary part of £ 

Magnitude of £ 

Derivative of ZI with respect to cp 
Derivative of ZR with respect to cp 
Real part of £ 

Value of ZI at the throat 
Value of ZR at the throat 
Value of ZPI at the throat 
Value of ZPR at the throat 
Square of the magnitude of £ 
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Table A~3« Input Parameters 


Variable 


Definition 


GAM 

CM 

SVN 


WC 

DWC 

NWC 

ANGLE 

RCT 

RCC 

IP 

AF 


Ratio of specific heats , Y 

Mach number at the nozzle entrance 

dJ (x) 

Nth root of the equation • — - = 0. Corresponds to 

Q.3C 

S . Values of S axe given in Table 1 for various 
acoustic modes 
Initial value of u> 

Increment of frequency 

Number of frequency points desired 

Nozzle half-angle, degrees 

Radius of curvature at the throat nondimens ionali zed 
with respect to the chamber radius 

Radius of curvature at the nozzle entrance nondimens ional- 
ized with respect to the chamber radius 

= 0: nozzle admittances are printed 

^ 0: pressure magnitude and phase are printed at each 

point along the nozzle 

Temporal damping coefficient \ 
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Table A-4. Output Parameters 


Variable 

Definition 

WC 

Nondimens ional frequency, u) 

YE 

Real part of the admittance as defined by Crocco in 
Ref. 2 

YI 

Imaginary part of the admittance as defined by Crocco 
in Ref. 2 

¥ 

Nondimens ional frequency 

SYR 

Real part of the specific admittance y 

SYI 

Imaginary part of the specific admittance y 
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Table A-5® Listing of the Computer Program Used to Determine 
the Ir rotational Hozzle Admittance (Page 1 of 10) 


1* 

2* 

3* 

4* 

5* 

6* 

7* 

8* 

9* 

10* 

11 * 

12 * 

13* 

14* 

15* 

16* 

17* 

18 * 

19* 

20* 

21 * 

22 * 

23* 

24* 

25 * 

26* 

27 * 

28 * 

29* 

30* 

.31* 

32* 

33* 

34* 

35* 

36* 

37* 

38 * 

39* 


common/xl/sam* svn, angle* rct* rcc /x2/t*rt* q* ri» r2* ip* wc.af 

C0MMCN/X3/Z1R* Z1I 
COMM0N/X4/ CM 
SAM = 1,233 

AF s 0 
IP=0 
RCC = 1 

RCT = 5. 437* 2/1 1,82 

N n C = 40 

Oi.C z 0*05 

ANGLE = 20 

CM = .25 

DO ICO I = 1*2 

I F ( I , E Q , 2 ) SO TO 5 

SVN s 0 
Nv.'C = 27 
GO TO 20 
5 SVN s 1.04129 
NWC = 20 
20 CONTINUE 

DO 2oO J = 1*3 
AF s 0,05*(J-2> 

IF (1 ,E0.2> SO TO 25 
WC = 0,55 
GO To 30 
25 WC = 1,55 
30 CONTINUE 

IF ( ip ,£Q, 0) SO TO 10 

WRITE <6* 1000) CM* SVN* SAM* ANGLE* RCT* RCC 
io call N0ZADM(CM* N W C* DWC ) . 

20o CONTINUE 
10o CONTINUE 

1000 FORMAT (46X* 28HPRESSURE MAGNITUDE AnD PHASE* //* 38X* 

1 14HMACH NUMBER = * F3,2* 7H SVN = * F6.4* 9H G&M'-IA = * p3 ’ 

2 . /* 22X* 15HN0ZZLE ANGLE = * F4.1* 2lH RADII OF CURVATUR 

3 * 9HTHR0AT = * F6.4* 12H ENTRANCE = » F6.4* //, 46x, 

4 2H X* 7X * 4HPMAG* 10X* 4HPARG* /) 

STOP 
END 


3^ 
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1* 
2 * 
3* 
4 * 
5* 
6 * 
7* 
6 * 
9* 
10 * 
11 * 
12 * 
13* 
14* 
15* 
16* 
17* 
16* 
19* 
20 * 
21 * 
22 * 
23* 
24* 
25* 
26* 
27* 
28* 
29* 
30* 
31* 
32* 
33* 
34* 
35* 
36* 
. 37* 
36* 
39* 
40* 
41* 
42* 
43* 
44* 
45* 
46* 
47* 
46* 
49* 
50* 
51* 
52* 
53* 
54* 
55* 
56* 
57* 
58* 
59* 
60 * 


Table A-5* Continued (Page 2 of 10) 


SUBROUTINE NOZADMtCMf NwC, DwC) 

DIMENSION DY (5»4) , 6(5) , GP<5> , Y(5) 

C.0MH0N/xl/3A-MrSVN.»ANSLE*.RCT#RCC/X'2/T».RT»:8»Rl»R2».IP»WC, AF 
COMMON// 3/ZlRfZlI 
DP = -0.001 

t = 3.1415927 * angle / i8o 

' *PITE(6f 1000) CM* SVN, 6AM, AF, ANGLE, P.CT , RCC 


00 1C NS 1, NWC 
20 ftC - WC + DWC 

25 RT =(c 1**0,5>*( (1* (GAM- 1 )*Cm*Cm/ 2 )**( <-6AM-1)/(4*(gAM-i) ) i 

1 > * ( (2/ (GAM+i ) ) ** { ( -SAM-1) /( 4* ( GAM-1) ) ) ) 

0 = (0«25*RT)+( (2/(GAM+l))**( (GAM*1)/(4*(GAM”1) ) ) ) 

PHIR = 1 

PHI I = 0 

R1 = RT + RC T * ( 1 - COS(T)) 

R2 = 1 - RCC*(i - COSCT)) 

R = RT 

P = 0 

U - 2 / (GAM+l) 

SRTR = ( RT * R C T ) * * 0 * 5 
AIR = -4 /( (6AM+l)*SRTR) 

BOR = -AIR +4*AF/(GAM+1) 

BO I = 4 * WC / ( gAM+1 ) 

SVNR s SVN/RT 

COR = WC * WC -{ ( SVNR*SVnR ) * 2 / (GAM+1) ) 

1 - AF*AF - 2* AF* ( GAM-1 ) / ( < GAM+1 ) *SRTR ) 

COI = -2 * WC * (GAM-1) / ( lGAM+l)*SRTR) - 2*AF*WC 

B1R = (24 ♦ 4*Gam)/I3*RCT*RT*(GAM+1) ) - 6*AF/(SrTR*<GAM+ 1 ) ) 

BIZ =6 * WC / (SRTR* (GAM+1 > ) 

C1R = 2 * (SAM _ 1) * SVNR * SVNR /(SRTR * <GAM+i)) 

1 - AF* (B1r+8*AF/(SRTR*(GAm+1) ) )* (GAM-1 ) *0,5 

Cl I = -31R * WC * (GAM - 1) * 0.5 
ZOR = (BOR*CoR + BOI+COl) / O0R*90R + BOI+BOI) 

ZOI = (B0R*C0I - BOI*COR) / <0OR*Bf)R + BOI*OOI > 

FI = B1R*Z0R - B.1I *Z01 - ZOR+ZOR+A 1 R + AlR*zOl*ZOl - ClR 
F2 = BlI*ZOR + B1R*Z0I - 2*AlR*Z0I*Z0R - Cli 
Z1R = ( F 1 * < A 1 R - BOR) - F2*B0I> / ( < A1R-B0R) * ( AjR-BOR) + 

1 B0I*B0I) 

Z1I = (F2* (AIR - BOR) ♦ F1*B0I) / ((A1R-B0R)*(Air-B0R) + 

1 301*301) 


C 

G ( 1 ) 
G(2) 
6(3) 

6 ( 4 ) 

G ( 5) 
OY ( 1 , 1 ) 
Oy (2, 1 ) 
DY ( 3, 1 ) 
DY( 4, 1 ) 
0Y( 5, 1 ) 
IQZ 
00 30 


= U 
= U 
= zOR 
- ZOI 

= PHIR * ZOr - PHI I * 201 
= PHI I * ZOR * ZOI * PHIR 
= -AIR 

= ZiR 
= ZH 
= PHIR 
= PH 1 1 
= 2 

I = 2,4 

call rxt?(5,op»p,g,sp,iqz) 

P = p + op 
U = G ( 1 > 

ZR = G(2) 

Zl = 6(3) 

PHIR s G(4) 

PHI I = G(5) 
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Table A-5. Continued (Page 3 of 10) 


61 ♦ 


0Y( 1 , 1.) = GPU) 

6a® 


DY(2,r) = GP (2) 

63* 


DY(3,I) = GPU) 

64* 


DY(4,D = GP<4) 

65* 

30 

D Y ( 5 , 1 ) s GP(5) 

66* 


Y(l) = U 

67* 


Y(2) = ZR 

68* 


Y(3) = ZI 

69* 


Y ( 4 ) = PHIR 

7C* 


Y < 5) = PM 1 1 

71® 


CALL ZAdAMS ( 5, OP, p, Y,DY, IOZ) 

72® 


IF (IP ,EO. 1) 30 TO 10 

73* 


U = Y<1) 

74* 


ZR = Y ( 2 ) 

75* 


ZI = Y ( 3 ) 

7b® 

P-'IR 

= ret) 

77* 

PmII 

- Y(‘j) 

78* 


GEAR = U* *0* 5 

79® 


C = 1 - U*C«5* (GAM-1) 

80® 


RMO = r.**( 1 /fGAv-l) ) 

81* 


F = OJAR / (GAM*RH0) 

62® 


Ic'IjZ ,r.O, 1) So TO 35 

83* 


ZJN = (U+ZR+AF ) * (U*ZR+Af) + <wC+U*ZI>*<WC+U*ZI> 

84® 


YR = -(ZR*(0*2r+AF) + ZI * < WC+U*Zl ) ) *F/ZDN 

85® 


Y I = " * ( w C * Z R - AF*ZI)/ZDN 

8b* 


GO T ) 4 0 

87 * 

35 

TR = Y ( 2) 

86 * 


T I = Y<3> 

89* 


TON = (|J + AF4TR-wC*TI)*(U+AF*TR-WC*TI) + <WC*TR)*(wC*TR) 

90* 


rR = -r * < U-aC*tI <- Afr *TR) /TON 

91® 


YI = r * ( .;C *TR + AF*TI ) / TON 

92* 


YI = f * m C * tr / TON 

93* 

4o 

SYR = r-AM* (C**< (.gAMfl ) /(2* (GAm-1) ))) *YR 

94* 


SYI = GAM*(C*» ( (gAM+1)/(2*<3AM-1) ) ) >*YI 

9b® 


W “* WM ^ :(C ♦ | 5) 

96* 

5C 

WRITE (6, 1005) WC, YR» YI, W* SyR* SYI 

97* 

lo CONTINUE 

96* 

100o F OR M / T 1 1 H 1 * u5X » 30HTHE0RETICAL NOZ..ZLE ADMITTANCES, //, 25X* 

99* 

1 

14RM.-1CH MUMPER = , F3..2, 7H SVN = , F6.4, 9H GaMMA s » p3,l 

100® 

1 

,2H decay coefficient = , F 6 . 4 , //, 

101® 

2 

22X* 15HN0ZZLE ANGLE = , F4,i* 2X. 21HRADI I OF CURVATURE? 

102® 

3 

, 9HTRR0AT = , F 6 * 4 , 12H ENTRANCE = * F6.4, //» 34X, 2HWC* 

103® 

4 

7x# 2-iYR* 8X* 2 l lYI • BX, 1 HW , 8X, 3HSYR* 8X, 3HSYI* /) 

104® 

1005 FORMAT ( 31X* F6 .4, 5F10.5) 

lOb® 

Rr.TjRN 

10 6® 

E'O 
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Table A-5® Continued (Page 4 of 10) 


I* 

2$ 

3* 

4* 

5* 

6 * 

7* 

6 * 

9* 

10 * 

11 * 

12 * 

13* 

14* 

15* 

16* 

17* 

18* 

19* 

20 * 

21 * 

22 * 

23* 

24* 

2b* 

25* 

27* 

28* 

29* 

30* 

31* 

32* 

33* 

34* 

3b* 

35* 

37* 

38* 

39* 

40* 

41* 

42* 

43* 

44* 

45* 


subroutine rktzinu, h, ti» u» oum» jOPT) 

C0MM0N/X2/T»RT»Q»R1»R2,IP*WC»AF 
DIMENSION U(5), A(5)> UZ(5)« FZ(*»»5) #DUM<5) 

A(l) = 0 
A t 2) = 0 
A<3> = 0.5 
A«4) = 0,5 
A<5) = 1,0 
TZ = T1 

DO lc J s 1, NU 
UZ(J) = u<u> 

10 DUM(J) = FZ(1,J) 

IF(J0PT ,EQ, 2) GO TO i5 
CALL RKTDIF( T Z,UZ»DUM) 

GO TO 20 

15 CALL RK2DIF(tZ,UZ»DUM) 

20 00 25 J : If NU 
25 FZ(l.J) = DUM(J) 

DMi I : 2.4 

TZ 5 Tl * A(I+i)*H 
DO 35 J = 1, NU 

UZ(j) = U ( J) + A(I+l)*H*FZ(J~lyJ> 

35 DuM(j) = F?(I , j) 

IF(J0PT .EQ, 2) SO TO 40 
CALL P<TDIr (TZ»Uz,DUM) 

30 TO 45 

*♦0 Call R<?DIF(Tz,Uz,DUM) 

<*5 DO 50 J = 1. NU 

5o FZ(I.J) = OUM(j) 

3o CONTINUE 

DO 55 J = 1, NU 

55 U(J) = u«U) + H*(FZ(1.J)+2*(FZ(2,J)+FZ(3»J) )+FZ(4*U>) / 6.0 

50 TO ( oO> 65 J » JOPT 
60 CALL R<TDIP(TZ.U »OUM) 

30 TO 70 

65 call ;«zdif<tz.u »dum> 

7 0 iF(IP.tO.O) 30 TO 75 

PR = y, C*U(*i)- Ull)*0UM(4) - AF*U(4> 

PI = -*C*u<4> - l.l(l ) *DUM(5) - AF*u(5) 

P'V .AS = SORT (ot<»PR + Ri*PI) 

PARS - AT A‘« { Pi/PR ) 

.v<IT- (ft* 1000) TZ, PM A 3 , PARS 
100o FORM; T(46X* cq . 4 > IX, F10.5» 3Xy Flo. 5) 

7=, Rc.TUAN 
EO 
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Table A-5* Continued (Page 5 of 10) 


it 

2 * 

3* 

4# 

5 * 

6* 

7 * 

at 

9 * 

lot 

lit 

12* 

13 * 

14 * 

15* 

16 * 

17 * 

18 * 

19 * 

20 * 

21 * 

22 * 

23* 

24* 

25* 

26* 

27* 

28* 

29* 

30* 

31* 

32* 

33 * 

34* 

35* 

3 b* 
37* 
38* 
35* 

4 1 ! * 
41* 
42* 

4 3* 


SUBROUTINE RKZDIF<P#G,gP> 

C0MMCN/Xl/GAv*,SVN#ANG|_E*RCT#RCC/X2/T»RT*a#Rl#R2#IP#wC,AF 

COMV10N/X3/Z1R.Z1I 

DIMENSION G ( 5 ) # GP(5) 

U s G { 1 ) 

ZR s G(2) 

ZI = 5(3) 

PHIr s G<4> 

P^II = G(5> 

IF(P) 15# 10, 15 

iO GP( 1 ) = 4/( (gAM+ 1 )*( <rcT*RT>**0.5> ) 

GP( 2 ) s ZiR 
GP(3) s Zll 
GP ( 4 ) = ZIR 
GP<5) r Zll 
GO To 20 

15 C = 1 - (GAM - 1) * U * 0,5 

R = 0 * ( (C)#!*(-» 1 /( 2 *(GAM» 1 ) ) ) ) ® (u#*-0«25) *4,0 
IF(R-l) 22# 22 # 50 
22 IF(R - Rl) 25, 30# 30 

25 DR = -< (2*rCT*<R-RT) - ( R-RT) + ( R-R'T) ) ** 0 ,5)/ (RT+RCT-R ) 

GO Tr 45 

30 IF ( R-R 2 ) 35# 40# 40 
35 DR = -TAN( T ) 

'30 TO 45 

4p DR = ((2*RrC*(l-R) - (R-l ) * ( R-l ) ) **0 , 5) / ( 1-R-RCC ) 

45 Dj = -(O**0.75)*(C** ( (2*GAM-.l)/(2*(GAM-l>>) )/<Q*(l-<GAM+l)*U*.5) 
1 ) 

SP( 1 ).= Du*DR 
3 ^ T C S 6 
5o GP( 1 > r 0 
5‘j A = J*(C-!|) 

BR = U* 3P ( 1 ) /C + 2*AF*U 
31 = 2 *.rtC*.| 

cr = #.:*ac » 5vn*sv;n*c/ (r*R) - af*af 
1 - ( 3 Av» i ) * AF*U* sp ( 1 ) *0,5* ( 1/C ) 

Cl = ) */,C*U*5p( 1) *0,5*( 1/C) - 2 *AF*wC 

5'(?)= ((:lR*ZR - BI *2 1 - CR) / A) - ZR*ZR + ZDZI 

G-\<3)= <(JI*ZR + 3R*ZI - CD / A) - 2*ZR*Zl 

3 ;J (4)= ZR*PhiR - ZI * a Hl I 

5P(b)= ZR*PHII + ZI*PhiR 

2q R : : TJ-N 
£ O 
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Table A-5* Continued (Page 6 of 10) 


1* 

2 * 

3* 

4* 

5* 

6 * 

7* 

B* 

9 * 

10 * 

11 * 

12 * 

13 * 

14 * 

15 * 

16 * 

17 * 

16 * 

19 * 

20 * 

21 * 

22 * 

23* 

24* 

25* 

‘26* 

27* 

26* 

29* 

30* 

31* 

32* 

33* 

34* 

3b* 

36* 


SUBROUTINE RKTDIFIP.G.gP) 

COMM0N/x1/GAM.SVN»ANGlE*RCT.RCC/X2/T*RT.G.R1.R2»IP,wC,AF 
DIMENSION G ( 5 ) » SP(5) 

U = G(l> 

TR s G ( 2 > 

TI = G(3) 

PHIR 5 G ( 4 ) 

PHI I B G(5) 

C = 1 - <gAM-1)*U*0,5 

R - Q * ( (C)**(-1/(2*(GAM-1) ) ) > * <U**-0.25) *4,0 
IF(R-l) 22.22.50 
IF(R-Rl) 25. 30. 30 

25 DR B -( (2*rCT*(R-RT) - (R-RT) * « R-rT) ) **0 .5) / < RT+rCT-R) 

GO TC 45 

30 IF(R-R2> 35.40.40 
35 DR = -TAN( T ) 

GO TO 45 

4o DR = ( (2*RcC*(l-R) - (R-l ) * (R-H > **0. 5) / < 1-R-RCC) 

45 D|j = -<U** 0 . 75 )*(C**( (2*GAM-l)/<2*< GAM- 1 ) ))> / (Q*( 1 -<GAM +1 ) *U* 


1 0.5)) 

GP(l)= 0U*DR 
GO TO 55 
5o GP(1) = 0 
5s A = U*(C-u) 

BR = U* DP ( l ) /C + 2*AF*U 
31 - 2 ♦ . , iO*u 

CR = wC*‘VC - SVN*SVm*C/(R*R> - AF*AF 
1 -( ) *AF*U*GP ( 1) *0,5* 1 1/C ) 

Cl = - (GAm- 1 )*wC*U*SP(l)*0,5*C 1/C) - 2 *AF*WC 
G?(2)= 1 - ( ,aR*TR-BI*Tl) - (CR*(TR*TR-TI*TI)- 2 *CI*TR*TI> )/ 
GP(3)= <-jR*tI - BI*TR + CI*(TR*TR-tI*TI) + 2*CR*TR*Tl) /A 
T2 = TR*TR * TI*TI 
S?<4)= (TR*P-lIR - Tl*PHlI)/T2 
.Pi 5>= ( TR* PRl I + TI*PHIR)/T2 


R-TU.N 


£f*D . 
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Table A-5. Continued (Page 7 of 10) 


!♦ 

3* 

§* 

6* 

7* 

a# 

9# 

10 * 

u# 

12* 

13# 

14* 
lb* 
16* 
17* 
1ft* 
19* 
20 * 
21 * 
22 * 
23* 
24* 
26* 
26* 
27* 
26* 
29* 
30* 
31* 
32* 
33* 
34* 
36* 
36* 
37* 
36* 
39* 
40 * 
41* 
42* 
43* 
44 * 
46* 
46* 
47* 
4e* 
49* 
50* 
51* 
52* 
53* 
54* 
55* 
56* 
57* 
56* 
59* 
6C* 
61 * 


SUBROUTINE ZADAM5(M*Hex*Y#DY#XQZ> 

C0MMCN/X1/G A*U SVN# ANGLE ♦«CT#RCC/X2/t*RT*Q*R1»R2* IP, wC,AF 

C0MM0N/X4/ CM 


DIMENSION COR ( 5) » DP(5>» DY<5»4)» PrEd< 5>» Y(5)» 8(5), 6P(5) 
lo CONTINUE 

DO 15 j - l»N 

PRED ( I ) = y(IJ+H*<55,*DY{If4)-59,*DY(I,3)+37,*DY(I»2J-g, + OY(i.i) 

1 1/24 eO 

15 CONTINUE 

X = X-H 
U = PRED(l) 

ZR = PRED (?) 


ZI = PRE0«3) 
PUR = PR El ( 4 ) 
PHI I = PRE3C5) 


C = 1 - (3AM-1)*U*o,5 

R = 0 ♦ ((C)**(-1/(2*(GAM-1)))) * (U**-0.25) *4,0 

iP(R-i) 1 7 * 1 7 » 1 00 
17 IP(R-Rl) 20, 25* 25 

2q OR = -< <2.rCT*<R-RT)-<R-RT)*<R“RT) )**0.5) / (RT+RCT-R) 

SO T' u C 

25 I-(R-R?) JO, 35, 35 
2 - '•» - »T AM ( 7 ) 

GO TO 40 

35 OR = ( (2*RT.C*< I-R) - (l-R)*Cl-R))**0.5) / U-R-RcO 
40 DU - -('J**0.75)*(C**( (2*GAu-i)/<2*<SAM-1))))/(Q*(1-<GAM+1)*U*0.5 
1 ) ) 

0P(1)= DR.Dj 

A = U* ( C~),J j 

3R = U*JP(1)/C+ 2*AF*U 
51 F 2 ♦ „C * 1 1 

CR = ftC*..C - ( SVN*SvN*C ) / ( R*R ) • AP*AF 
1 * ( 3 Av_l ) * AP*U*DP ( 1 ) *0, 5/C 

Cl = - ( 3 A' 1 -! ) *riC*U*Qp ( 1 ) *0* 5/C - 2*AF*WC 
DP<2)= ({jR./.R - BI*ZI - C R ) / A ) - ZR*ZR ♦ 2 l*Zl 
D 3 <3!= (()I*?R + BR*ZI - CI)/A) - 2*ZR*ZI 
OP(4)F ZR. 3 HJR ~ ZI *P-tI I 
0P(5)= Z R ♦ P H 1 1 * ZI *P-UR 
DO 45 I = 1,N 

Cc»R ( I ) = Y(I)+H*(0y ( I,2)-5,*0YlI,3)+19,*DY<I»4)+9.*0P(I))^24.0 
45 Ytl) = ( 251 , +C0R ( 1 ) ♦ 19« *PRED ( I ) ) / 270. 

U = Y( 1 ) 

ZR = Y<2> 

ZI = Y ( 3 ) 

PUR = Y ( 4 ) 

P-U I = Y ( b ) 

C = 1 - (3AM-l)*U*o.5 
52 DO 5fi I = l*N 

DY ( I , 1 ) = DY { I , 2 ) 

Dr ( 1 ,2) = D Y ( I * 3) 

55 DY (1,3) = D Y ( I , 4 ) 

Z'-'63 = ( ZR*ZR + ZI*ZI)**0.5 
IP ( V A 3 > 10 ) 60 * 90, 90 

60 R = 0 * ( (C)**(-1/(2*(GAM-1)))> * (U**-0.25) *4,0 
IP.(R-I) 62» 62, 100 
62 IF(R-Rl) 65,70,70 

65 DR = -( (2*rCT*(R-RT) - (R-RT)*<R-rT))**0.5)/(RT+RCT-R) 

SO T- 65 

7c t c (R-R2) 75,80,80 
75 DR = - T A N ( r ) 

3.0 Tu 65 
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Table A~5® Continued (Page 8 of 10) 


62* 

63* 

64 * 

6b* 

66 * 

67* 

63* 

69* 

7U* 

71* 

72* 

73* 

74* 

75* 

76* 

77* 

78* 

79* 

80* 

81* 

82* 

83* 

84* 

85* 

86 * 

87* 

88 * 

89* 

90* 

91* 

92* 

93* 

94* 

95* 

96* 

97* 

98* 

99* 


®0 

®5 


67 

90 


y-i = ( (2*?nC*(l-R) - <1-R)*(1-R> )**0,5)/{l-R-RCC) 

DJ = -(U* *0.7.5) *<C** ( <2*SAM-1)/(2*(GAM-1> ) )>/(Q*<l* 
DY(1,4)= 0R*hU 
A = U*lC-u) 

BR = u*'DV(i*4)/C + 2 * AF*U 
Bl s 2*#C*y 

CR = WC*WC - C 5 VN * S VM * C ) / ( R * R ) - AF*AF 
l -(SAM-1)*AF*U*Dy(1»4)*0,5/C 

Cl = -<SAM-1 )*wC*U*DY<1»4)*0,5/G -2*aF*wC 
07(2.4)5 (BR*ZR - 9 1 *2 1 -CR)/A - ZR*ZR + ZI*ZI 
DY(3,4)r (3l*ZR + BR*ZI -CD/A - 2*ZR*Zl 
DY (4*4) s ZR*pHIR - ZDPHII 
DY(5,4)= ZR*PHII + ZI+PHIR 
IF (IP ,EQ. 0) GO TO 87 

PR s WC*PHH - U*DY(4*4) - Ap*PHlR 
PI s -wC*PhIR -U*DY (5*4) - AF*PhH 

PM A 5 = (PR*PR + PI*PI)**,5 
PARS = A TAN (pi /PR') 

WRITE (6* 1000) X* PMAS, PARS 
GO TO 10 
IQZ = 1 

Z2 = ZMAG+zMAS 

Y<2) = ZR/Z2 
Y (3) = -ZI/ZR 
ZPR = OY 12*4) 

ZPI = 0Y(3*4) 

OY( 2 > 4) s -<ZPR*(ZR*ZR 


DY(.3,4)s (2*zPR*ZR*ZI 


G(l> 

S<2) 

G<3> 

S.C 4) 

0(5) 
DYU, 
DY(2,i)= 
D Y (3*1 ) s 
OY (4*1)5 
DY ( 5 , 1 ) s 


U 

= Y (2) 

= Y ( 3) 

= PHIR 
= PHI I 
1)= D Y ( l * 4 ) 
DY( 2*4) 
OY ( 3 * 4) 
PHIr*ZR 
PHII+ZR 


„ ZI*ZI) ♦ 2*ZR*ZI*ZPI)/(Z2*Z2> 
- ZPI* (ZR*ZR - ZI*ZI))/<Z2*Z2) 


- PHI I* Z I 
+ PHir*ZI 


100* 

00 95 I = 

2*4 

101* 

call 

RKTZ<5* 

102* 

X 

s 

x+h 

103* 

u 


sd» 

104* 

TR 

- 

G(2> 

105* 

TI 

s 

■0(3) 

lOfa* 

PHIR 


S ( 4) 

107* 

phi i 

- 

S ( 5) 

108* 

DY ( 1 * I ) 

s 

GP(1) 

109* 

DY ( 2 * I ) 

" 

GP(2) 

110* 

DY ( 3* I ) 

s 

GP<3> 

111* 

DY(4* I ) 

s 

GP(4) 

112* 

95 Dy ( 5* I ) 

” 

SP(5) 

113* 

Y(l) = U 



114* 

Y<2) = TR 



115* 

Y(3) - TI 




Ufa* Y(4) = PHIR 

117* Y<5) = PHI I 

118* CALL TADAMS(n*H*X*Y*Oy*IQZ*IQ) 

119* GO TO <10* loO)*IQ 

120* 100o FORMaT(46X*F6.4»1X,F10.5*3X»F10.5) 

121* 10o RETURN 

122* END 


(GAM+1) *U/2) ) 


41 



Table A-5® Continued (Page 9 of 10) 


1* 

3* 

4$ 

5# 

6* 

7* 

8 * 

9* 

10 * 

11 * 

12 * 
13* 
14* 
15* 
16* 
17* 
16* 
19* 
20 * 
21 * 
22 * 
23* 
24* 
2b* 
2 b* 
27* 
2b* 
29* 
30 * 
31* 
32* 
33* 
34* 
3b* 
36 * 
37* 
38* 
39* 
4Q* 
41* 
4 2 * 
43* 
44* 
4b* 
46* 
47* 
4 6* 
4?* 
50* 
51* 
52* 
53* 
54* 
5b* 
56* 
57* 
58* 
59* 
60* 
61* 
62* 


SUBROUTINE TADAMS(NfH P x»Y#Oy#lQ2»lO) 

CDmwn/X1/GAM,SVN# angle# RCT ? RCc/X2/t»RT^0*R1jR 2» IP, !afC,AF 
C0MM0N/X4/ CM 

DIMENSION COR (5) e DP<5) * 0Y<5r4>' PRED (5) * Y (5) * S<5>, GP<5> 
lo CONTINUE 

00 15 I = 1 » N 

PREO(I) = Y<I>+H*(55*0y(I»4)-59.*DY(I»3)+37,*0Y(l*2)-9*DY(I»l>)/ 
i 24,0 

i5 CONTINUE 

X = X+H 
U = PRED(l) 

TR = PREO ( 2 ) 

T I = PRE0<3) 

PH I R = PREDU) 

PH I I = PREDtS) 

C = 1 - (-,AM-l)*U*.5 

R = 0 * * (U**-0.25) *4,0 

IP(R-l) 17-17*100 
1 7 IEn-71) 20, 25 1 25 

20 Dr = -( (2*rCT*(R-RT) - (R-RT) * (R-rT ))**.. 5)/ (RT + RCT-R) 

50 T 40 

25 IF { R-R2 > 30, 35, 35 
3o 3'< = -TAN(x) 

3D .T c 4 0 

35 DR - ( (2* R'C*(1-R> - (1-R)*M-R))**.5)/(1~R-RCC) 

4 0 Dj =-( J«*,?5)* (C** ( { 2+GA'l-l ) / { 2* (G^M-1 ) ) ) )/(0*U-(6 aM<-1) *U*.5> ) 

D :> (1) = DR*Dj 
A - u*(C-U) 

3r = U>J a < 1 ) /C+ 2*AF*U 
5 I = 2 * a C * • ) 

CR = *C»*D ” ( SVN*S'7N*C ) / ( R*R ) - AF*ftF 
1 - ( 3* “—1 ) * AF*U*DP( 1 ) *0,5/C 

Cl r - ( 5A>'«1 ) *wC*U*Dp ( 1 ) *0 . 5/C - 2*Ap*WC 
D p ( 2 ) = 1 4 ( - •RiTri + DI*Tl + Cf l*(TR*TR~Tl*TI )-2*CI*TR*Tl)/A 
0^(51= ( - .1 R ♦ T I - B I * T R + CI*(TR*TR - TI*TI ) + 2*CR*TR*TI ) /A 
T 2 = T R ♦ T R + TI*TI 
DP ( 4 ) = <TR* p> lIR - Tl*?HlI >/T2 
DP ( 5 ) = ( T R * P -1 1 1 + TI*PHIR>/T2 

JO 4 -. I = 1#N 

Cc R(I) = Y<I) + H*«DY(I*2)-5**DY<i,3)+19.*DY<I»4)+9,*DP(I))/2»1.0 
45 Y(I) = (251, *C0R( I ) + 19. *PRfe.D{ I ) ) /270, 

U = V ( 1 > 

TR = Y l 2 ) 

II - Y(3> 

PUR = Y ( 4 ) 

P 'II = Y ( 5 ) 

C = 1 - t GAM-1 ) *U*,5 
S 2 DO S 1 I = 1,‘i 

D Y 1 1 , 1 ) = D Y ( 1 , 2 ) 

D Y ( I » 2 ) = OY ( I f 3) 

55 DY< I, 3) = D Y ( I # 4 ) 

T 2 = TR«TR +■ TI*TI 
T''Aj = T 2 * * • 5 
IF(T Ac, - 10 ) 60- 90, 90 

60 R = Ci * ( ,C) + +I-1 /{ 2’»(SAM-1)) ) ) * (u**-0,25) *4.0 
IF { r- 1 ) 62- 62, 100 
&2 I-(R-Rl) o5, 70,70 

65 Dr = - ( (2*rCT*(R-RT)-(R~RT)*(R-RT) >**.5)/<RT+r0T-R) 

50 Tc 65 

7o I~( R-R2 ) 75,rO,0O 
?5 DR = -TAN(t) 

50 T- 35 
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Table A-5® Continued (Page 10 of 10) 


63$ 
64 * 
66$ 
66$ 
67$ 
66$ 
69* 
70* 
71* 
72* 
73* 
74* 
75* 
76* 
77* 
76* 
79* 
80* 
81* 
82* 
83* 
84* 
85* 
86 * 
87* 
88 * 
89* 
90* 
91* 
92* 
93* 
94* 
95* 
96* 
97* 
98* 
99* 
100 * 
101 * 
102 * 
103* 
104* 
105* 
106* 
107* 
108* 
109* 
110 * 
111 * 
112 * 
113* 
114* 
115* 
116* 
117* 
118* 
119* 
120 * 
121 * 
122 * 
123* 


8 G D-< = ( (2*^C*<1-R) - <1-R)*<1-R> )**«5)/(l-R~HCC) 

85 dj = -(0**,75)*<C**( (2*SAM-i)/<2*(GAM-l))))/(Q*(l-CsAM+l)*U*,5) 
DY< 1,4)= DRoJ 
A = U*(C-ij) 

3* = UOY ( l » 4) /C + 2* AF*U 
3l - 2*wC*u 


CR = WC*WC - (SVN*SVM*C)/(R*R) - AF*aF 


1 -(GAm-1)*AF*U*DY<1»4)*0,5/C 

Cl = »(GAM-1 )*wC*U*DY( 1»4)*0,5/C „2*AF*wC 
DY<2,4)= 1 * <-BR*TR + BI*TI + CR*(TR*TR - TI*TI> - 2*CI*TR*TI ) /A 
DY<3,4>= (-BR*TI - BI*tR + CI*(TR*TR - TI*TI> ♦ 2*CR*TR*TI ) /A 
DY (4,4)= <TR*PH1R - PHII*TI)/T 2 
0Y(5,4)= (TR*PHII + PhIR*TI)/T 2 
If C IP , EQ. 0) GO TO 87 

PR 5 WCoPHII - U*DY ( 4 » 4) - AF*PHIr 

pi = -wc*phir -u*dy<5>4> - af*phii 

PMAS = <PR*PR + PI*PI)**,5 
PAR3 = ATAN(pI/PR) 

WRITE (6, 1000) X, PMAG, PARG 
87 GO TO 10 
9o IOZ = 2 
Y<2> = TR/T2 
Y<3> = -TI/T2 
TPR = DY(2»4) 

TPI = DY(3,4> 

0Y(2,4)= -<TpR*(TR*TR - TI*TI ) + 2 *tR*tI*TPI)/(T 2*T2> 
DT(3,4)=(2*TpR*TR*Tl-TPl*<TR*TR-TI*Tl))/(T2*T2) 

G<1) = U 
G<2) = Y ( 2 > 

G<3> = Y ( 3) 

G(4) = PHIR 
G<5) s PHI I 
0Y(1,1)= OY ( 1 , 4) 

DY(2,1)= i)Y(?,4) 

DY ( 3, 1 ) = OY ( 3 , 4) 

0Y(4,1)= (PHIR*TR - PHII*TI)/T2 
0Y(5,1)= (PHII*TR - PHIR*TI)/T2 


DO 95 I = 2,4 

call Rktz(5,h,x,g,gp,iqz) 

X = X+H 

U = 5(1) 

ZR = G<2> 

ZI = G 1 3) 

PHIR = G<4) 

phi i = 5(5) 

DY ( 1 , 1 ) = GP(1) 

0Y( 2, 1 ) = gP(2) 

DY (3, I) S S P(3) 

D Y < 4 1 1 1 = gP 14) 

95 Dy t 5,1) = GP(5) 


10o 

1000 

105 


Y<1) 

Y<2) 

Y<3» 

Y < 4 ) 
Y< 5) 
10 

GO TO 


= U 
= ZR 
= ZI 
" PHIR 
= PHII 
= 1 
105 


IQ = 2 
FORMAT (46X* 
RETJRN 


F6.4» IX, FlO.5, 3X, Flo. 5) 


END 
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Table A- 6.- Sample Output 


id to to co o N cr h c* or d m t- cv o- to p- 

o cv d d f**- co fo sC <x) *— t p- cr «—t tn or o fs* n? m 

‘Z 1 CM CC \D N vO D N W CM C OD tO W -f cp C 

ry vc *~mp cr ?o c* ovr- in d to -+• ve O' /v A * 

fOD (\iC\l HOC O C;r! (VIO^ IT; vD p' &' c. ^ 


f ? t 


I i l 


o ^ Ch n vO vO ro to ^ c cc o d n h c c ^ in 

r- in d m x ccogn lo d -* Xi cr o cp er o 

vo -W 0- d <h o crj *.0 fO p to x cr r- c\j -b — « cv cr 

K) M O P ro VO IT) ^ fO W C ^ m CVl ^ ^V1 vD P VC P 

rO lO K> <V OJ : CV CM <M CV CV] c\j ^ ^ ^ 0 0 o w jr> S) 


w o .co vO d cv o cr* p' lo to r-t cr *o \c zr eg *-• x r- 

co o *h m to h*- cr o cv d vo ro cr «h to ld rr A ^v 

m \0 vo vo vo vc sT N r- r- n p. or,« <7 a* cc .or a n. 

c-’ c in c: in o ir. c- in c m in o <n * m tn 

vO vC.! P- P* Cl! CO O'' 0 s - C ’ .O' r~' *H CV CV m rO ~i /-» I:"- 


h h (\i eg c\i cv eg cv cy oj pi tv cj cv 


?o m to tv r~i a* x n- vo cc ^ to »h ^hi rr> r-« cv 

co O' m d- oj \c d o th lo rg to <r d r. in p- 

u d •:.: in cr *h oj o d a, -g oj o n- cr p- in 

in «h r- eg r- to co n d cr, g ^ co vr m lo vo n cn 

fo mo* oj r - o o o *~! cv Cv to tr ^ ^ ^ 


i l l J 1 ff i i 


k; w c- if) cr- *-• a a, cv r- 0‘ r- m a> k - 

r^ o cv ^ \0 vO in O' m X) lo cr S n ^ 

oj o cd n- x» xo m d nr d n \0 cr- r- • ^ 
d: MTdfC WHOO to f- ID r; o s o - ir 

0 oj cm ou o’ Oj eg Cl t— T— »■-: »-w ^ r 

****♦«»*«•*•*•*** 

1 i » I I * 1 I t ! I 1 * ! I J 


o o o c.’ cr c.' 
c ■ c r: c? o c > 
■o til ' r- in o 5 J1 
\f> \0 ■ Nfc * n ;o 


C > C> C* C*‘ c 
p- to o in o 
:i' # T; x w 


r* c r* 
r n in 

H :g af 


vV ,V Ov; 0 i\l 
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SUBROUTINE 

NOZADM 



- 

STEPSIZE 
DP = -0.001 

, 


WRI 
CM, SVT 
RCC , AT 

:te 

r, RCT, 
fGLE 


) 


DO THROUGH 

a 

FOR N = l,NPTy 


WC = WC + DWC 
COMPUTE 
F, RT, Rl, R2 
WC 


COMPUTE INITIAL 
VALUES OF 
PHIR, PHII, U, 
A,B,C ,ClR,ClI, 
BlR,BlI,ZR,ZI 


ZR = 

ZOR 

ZI a 

ZOI I 

G(l) 

= 

u 

G(2) 

== 

ZR 

G(3) 

= 

ZI 

0(4) 

= 

PHIR 

0(5) 

SS 

PHII 

DY(1 

,1) 

to 

«»y«C!9 

HS1 

,JQ£L— 


3 


DO THROUGH 

P 

FOR I = 2,4 


I 


CALL 

RKTZ 


P = P + DP 
U = G(l) 

ZR = G(2) 

ZI = G(3) 
PHIR = G(4) 
PHII = G(5) 
COMPUTE 
DERIVATIVES 


COMPUTE 
NEW VALUES OF 
C, U, CM 


DY(1,I) 
DY ( 2 , 1) 

DY(3,I) 
DY ( 4 , i) 


= GP(1) 
= GP(2) 
= GP(3) 
= gp(4) 
S GP(5l 



COMPUTE 
QBAR, YR, YI 
FROM ZR, ZI 


COMPUTE 
QBAR, YR, YI 
FROM TR, TI 


WRITE 

WC, YR, YI, 
W, SYR, SYI 



RETURN 
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Figure A-l. Continued (Page 5 of 10) 





Figure A-l. Continued (Page 6 of 10) 
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COMPUTE DR FOR 
SECTION III 


COMPUTE DR FOR 
SECTION II 


COMPUTE DR FORLJ 
SECTION I P 


* 


COMPUTE 


C,BU,A,BR,BI, CR, 
Cl jPHIR jPHII , 
DY(1,4) THROUGH 



r 
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